% Computes bi-variate transfer functions between two channels, given
% two input & one output channel
% Latest date 21.11.2002
% bug fixed - correct signs for the cross spectra
% compared with readedi1 erros, z coh found correct

%zij = za zb - zc zd;
%tf = zxy zyx zxx zyy
function[tf] = TransFn(SPM1,ProcDef),

nfreq=ProcDef.nfrq;
EQN=ProcDef.equations;
neq=ProcDef.neq;

if isstruct(SPM1) == 1,
   for i = 1:nfreq,
	   SPM(i,:,:) = SPM1.spectra(i).data;
   end;
end;

tf = zeros([nfreq neq+2]);

for l = 1:neq,
i = EQN(l,1);
j = EQN(l,2);
k = EQN(l,3);
p = sqrt(-1);

za = SPM(:,j,j);
if i>k,
   zb = SPM(:,i,k)+p*SPM(:,k,i);
elseif k>i,
   zb = SPM(:,k,i)-p*SPM(:,i,k);
end;

if k>j,
   zc = SPM(:,k,j)-p*SPM(:,j,k);
elseif k<j,
   zc = SPM(:,j,k)+p*SPM(:,k,j);
end;

if i>j,
   zd = SPM(:,i,j)+p*SPM(:,j,i);
elseif i<j,
   zd = SPM(:,j,i)-p*SPM(:,i,j);
end;

den = SPM(:,j,j).*SPM(:,k,k)- (SPM(:,k,j) - p*SPM(:,j,k)).*(SPM(:,k,j) + p*SPM(:,j,k));
top = [za.*zb-zc.*zd SPM(:,k,k).*zd-conj(zc).*zb];

tf(:,l) = (top(:,1)./den);
tf(:,l+2) = (top(:,2)./den);

end;

for i = 1:nfreq,
   data(:,:) = SPM(i,:,:);
   c1 = cohe(tf(i,3),tf(i,1),data,'Ex');
   c2 = cohe(tf(i,2),tf(i,4),data,'Ey');
   dof = SPM1.spectra(i).dof.avgt*SPM1.spectra(i).dof.avgf;
   [tf(i,7),tf(i,5)] = BivError(c1,dof,data,'Ex');
   [tf(i,6),tf(i,8)] = BivError(c2,dof,data,'Ey');
end;

